The present invention generally relates to seismic exploration, and more particularly to accurate prediction of 3-D acoustic reverberations or multiples for the purpose of coherent noise suppression and improved interpretation of 3-D seismic data.
Seismic exploration involves the study of underground formations and structures. In seismic exploration, one or more sources of seismic energy emit waves into a subsurface region of interest, such as a geologic formation. These waves enter the formation and may be scattered, e.g., by reflection or refraction. One or more receivers sample or measure the reflected waves, and the resultant data are recorded. The recorded samples may be referred to as seismic data or a xe2x80x9cseismic tracexe2x80x9d. The seismic data contain information regarding the geological structure and properties of the region being explored. The seismic data may be analyzed to extract details of the structure and properties of the region of the earth being explored.
In general, the purpose of seismic exploration is to map or image a portion of the subsurface of the earth (a formation) by transmitting energy down into the ground and recording the xe2x80x9creflectionsxe2x80x9d or xe2x80x9cechoesxe2x80x9d that return from the rock layers below. The energy transmitted into the formation is typically sound energy. The downward-propagating sound energy may originate from various sources, such as explosions or seismic vibrators on land, or air guns in marine environments. Seismic exploration typically uses one or more sources and typically a large number of sensors or detectors. The sensors that may be used to detect the returning seismic energy are usually geophones (land surveys) or hydrophones (marine surveys).
During seismic exploration (also called a seismic survey), the energy source may be positioned at one or more locations near the surface of the earth above a geologic structure or formation of interest, referred to as shotpoints. Each time the source is activated, the source generates a seismic signal that travels downward through the earth and is at least partially reflected. Seismic signals are partially reflected from discontinuities of various types in the subsurface, including reflections from xe2x80x9crock layerxe2x80x9d boundaries. In general, a partial reflection of seismic signals may occur each time there is a change in the elastic properties of the subsurface materials. Reflected seismic signals are transmitted back to the surface of the earth, where they are recorded as a function of traveltime. Reflected seismic signals are typically recorded at a number of locations on the surface. The returning signals are digitized and recorded as a function of time (amplitude vs. time).
FIG. 1A illustrates a seismic source S at the earth""s surface, referred to as the free surface, generating reflected seismic signals from a sub-surface which are measured at four receivers, R1-R4, as shown. One should note that all free surface reflected energy acts as a secondary source of seismic signals from the points of reflection. FIG. 1B illustrates various examples of primary reflection raypaths, where the signal refracts at each boundary between layers, i.e., at each sub-surface reflector. Note that primary reflections do not involve downward reflections from the free surface or any of the sub-surface reflectors.
In seismic analysis, the term xe2x80x9cmultiplesxe2x80x9d refers to multiply-reflected seismic energy, or any event in seismic data that has incurred more than one reflection in its travel path. Depending on their time delay from the primary events with which they are associated, multiples are commonly characterized as short-path, implying that they interfere with the primary reflection, or long-path, where they appear as separate events. Multiples from the water bottom (the interface of the base of water and the rock or sediment beneath it) and the air-water interface are common in marine seismic data. The presence of multiples may obscure or interfere with primary reflection signals and may thus act as xe2x80x9cnoisexe2x80x9d when analyzing seismic data.
FIGS. 2A-2D illustrate a variety of multiple events, as are well known in the art. FIG. 2A illustrates water layer reverberations, where the signals are xe2x80x9ctrappedxe2x80x9d in the water layer between two strong reflectors, specifically the free surface and the bottom of the water layer. These multiples tend to decay slowly and obscure the primary reflection energy from deeper reflectors or sub-surfaces. FIG. 2B illustrates slightly weaker multiple events, referred to as xe2x80x9cpeg-legxe2x80x9d multiple events, characterized by an additional roundtrip through the water layer just after emission (source-side) or just before detection (receiver-side). FIG. 2C illustrates a variety of what are known as xe2x80x9cremainingxe2x80x9d surface-related multiple events, where the first and last upward reflections are below the first (water) layer, and there is at least one reflection at the free surface in between. These multiple events are typically weaker than the water layer reverberations and the peg-leg multiples, but may be considerable if a highly reflective structure (e.g., salt or basalt) is involved. Finally, FIG. 2D illustrates xe2x80x9cinternalxe2x80x9d or xe2x80x9cinterbedxe2x80x9d multiple reflections, which are generally much weaker than the surface-related multiple reflections of FIGS. 2A-2C. Internal multiple reflections have a downward reflection at a boundary in the subsurface. If the primary reflection events have reflection amplitudes on the order of the reflection coefficient r, then the first order surface-related multiples have amplitudes on the order of r2, and the first order internal multiples of r3, where |r| less than 1. However, in many common geologic settings of commercial interest in hydrocarbon exploration, multiple events will interfere with less energetic primary events from deeper reflectors, and even internal multiples may obscure primary reflections from deeper reflecting boundaries.
FIG. 3 illustrates the time relationship between the arrival of a source sample at a receiver R, referred to as the sample time, and the arrival times of a multiple reflection from some representative subsurface reflector. Note that time increases downward along the vertical axis shown. As FIG. 3 illustrates, the original sample produces the multiple event shown to the right of FIG. 3. These arrival times are those of a primary reflection from that subsurface reflector due to a source located at the receiver location, but delayed by the sample""s arrival time. Note that the sample may itself be a measurement of a multiply-reflected signal from a source located elsewhere, and that any energy reflected from the surface is, according to the principle of superposition, considered to be another source. In discrete terms, FIG. 3 illustrates the concept that the response of any sample of recorded upcoming energy is a delayed copy of a record from an impulsive seismic source signal at that recording location, appropriately scaled by the sample""s amplitude and the free surface reflection coefficient.
Multiple source activation/recording combinations may be collected to create a near continuous profile of the subsurface. In a two-dimensional (2-D) seismic survey, the recording locations are generally laid out along a single line. In a three dimensional (3-D) survey the recording locations are a really distributed across the surface. In simplest terms, a 2-D seismic line can be thought of as providing a vertically-oriented cross sectional picture of the earth layers as they exist beneath the recording locations. A 3-D survey produces a data xe2x80x9ccubexe2x80x9d or volume that is, at least conceptually, a 3-D picture of the subsurface that lies beneath the survey area. In reality, both 2-D and 3-D surveys interrogate some volume of earth lying beneath the area covered by the survey.
After the seismic data have been collected, the seismic data may be xe2x80x9cimagedxe2x80x9d, analyzed, or otherwise processed to produce a seismic profile or pattern indicating various characteristic structures or signatures which may correlate with known geological formations. Seismic processing generally refers to alteration of seismic data to suppress noise, enhance signal and migrate seismic events to the appropriate location in space. Processing steps typically include analysis of velocities and frequencies, static corrections, deconvolution, normal moveout, dip moveout, stacking, and migration. The migration step can be performed before or after stacking. Seismic processing facilitates better interpretation because subsurface structures and reflection geometries are more apparent. A number of seismic exploration techniques may be performed to analyze the seismic data, including pre-stack imaging, pre-stack migration and normal moveout.
Moveout refers to the difference in the arrival times or traveltimes of a reflected wave measured by receivers at two different offset locations. Said another way, moveout refers to arrival times of events with increasing source-to-receiver offset. Normal moveout (NMO) is moveout caused by the separation between a source and a receiver in the case of a flat reflector under a homogeneous overburden. NMO may also be considered an approximation of moveout by hyperbolic trajectories. Dip moveout (DMO) occurs as an effect in addition to NMO when reflectors dip. Inverse moveout refers to the process of predicting and restoring arrival times with offset. Inverse normal moveout is a simplified process that uses hyperbolas to restore the arrival times.
Convolution is a mathematical operation on two functions that represents the process of linear filtering. Convolution can be applied to any two functions of time or space (or other variables) to yield a third function, the output of the convolution. Although the mathematical definition is symmetric with respect to the two input functions, it is common in signal processing to say that one of the functions is a filter acting on the other function. The response of many physical systems can be represented mathematically by a convolution. For example, a convolution is commonly used to model the filtering of seismic energy by the various rock layers in the Earth.
As mentioned above, the presence of multiples may complicate the analysis of seismic data sets by creating reflection arrivals masking or masquerading as primary reflections for a formation of interest. Thus, various techniques have been developed for removing or at least ameliorating the effects of multiples from seismic data. An approach called surface-related multiple attenuation (SRMA), devised in 1-D at Stanford University in the late- 1970""s and extended and developed at Delft University of Technology in the mid-1980s, referred to herein as the Delft approach, is a prediction and subtraction process, where multiples are first predicted, then subtracted from the seismic data. The Delft surface-related multiple attenuation algorithm and its various extensions to subsurface interbed reverberations are considered by practitioners of the art of seismic data processing to be powerful tools for the prediction and attenuation of multiple reflections in 2-D prestack seismic data. The Delft surface-related multiple approach may be understood by means of an algebraic relationship. In the Delft approach, each shot record approximates the response of the subsurface to an ideal impulsive source at the shot location. Each trace (usually multiplied by xe2x88x921, the idealized free surface reflection coefficient) approximates energy being reflected downward from the free surface. Convolving the trace with a shot record whose source lies at, or nearby, the receiver location therefore approximates the response of the subsurface to the free surface reflected trace energy. This relationship may be expressed as follows:
[S+(xe2x88x92D)]*R=Dxe2x80x83xe2x80x83(1)
where S are the actual source signals, D are the upcoming recorded data, xe2x88x92D are the free surface reflected data, and R is the response of the subsurface to ideal impulsive sources without the free surface reflector. In other words, the recorded data are the response of the subsurface to ideal impulsive sources without the free surface reflector to the initial source illumination plus the downward-reflected data. Assuming one can measure or estimate the source field, S, one may solve this equation for R and thereby eliminate all free-surface multiples.
To highlight how powerful this concept is, FIG. 4 illustrates some of the complexities that may result from non-parallel boundaries between subsurface layers, where the primary reflection path from the source S0 to the receiver R1 is denoted by the heavy line labeled Recorded Data. Note the various other signal paths that propagate through reflection and refraction at the boundaries of subsurface 1 and subsurface 2. Note also that the signal reflected from the free surface at the location of the receiver R1 is considered to be a secondary source S1. This signal may be convolved with a shot coincident with the receiver location. In reference to FIG. 4, the Delft algorithm may be stated thus: each recorded trace is convolved with a shot record, possibly interpolated, from a shot at that trace""s receiver location. Superimposing these individual records across the whole survey line produces a prediction of the surface-related multiples for the line. These predicted multiples are then adaptively subtracted from the original data to suppress surface-related multiples.
It should be noted that most 2-D implementations of the Delft approach take advantage of the near-continuous linear recording geometry of normal 2-D seismic data acquisition by invoking the xe2x80x9cprinciple of reciprocityxe2x80x9d wherein a source activated at location S and recorded at receiver location R produces the same signal as the same source activated instead at location R and recorded by the same receiver repositioned at location S. In other words, sources may be mathematically interchanged with receivers, say in a computer, to effectively enlarge the recording geometry of the original near-continuous profile.
The relationship (1) holds true for 3-D as well as 2-D or 1-D earth models. However, in 3-D, recording cables need to be replaced by large areal recording arrays. In other words, to directly deploy the Delft approach in 3-D, data would have to be recorded over an areal grid (a 2-D grid, such as a rectangular grid) from shots distributed over the grid, i.e., shotpoints and receivers would be required over the entire grid. This type of field acquisition, while technically possible, is simply too expensive with current acquisition technologies to justify for everyday petroleum exploration. Thus, while in theory the Delft approach can be directly applied in three dimensions, the practical cost of acquiring a dense coverage of areal shot records has historically made this extension almost purely of theoretical interest.
Various researchers have attempted to adapt the Delft approach to swaths of parallel towed streamers where there is at least some limited areal receiver coverage, for example, by estimating some local crossdip and warping the data to compensate for the out-of-plane effects, by extrapolating the data to a wider and denser crossline coverage, or by implementing an anti-aliased summation within the towed streamer swath to compute a gracefully-incomplete 3-D Delft response. However, none of these attempts truly overcomes the severe aperture limits of real world data acquisition, and none has truly addressed the issues related to diffracted multiples, which are inherently of a fully 3-D nature.
Thus, improved systems and methods are desired for performing 3-D seismic analysis and processing of 3-D multiples.
Various embodiments of systems and methods for predicting or estimating 3-D multiples from 3-D stacked seismic data are presented. In one embodiment, prestack seismic data may be received, e.g., directly from a plurality of receivers, e.g., via wired or wireless network, satellite, telephone, delivered storage medium, or any other means of transfer, or from an intermediate system, e.g., a third party. A stack of data, i.e., stacked seismic data, and a corresponding stacking velocity field may be generated from the prestack data, where the stacked data includes a plurality of stack traces. Each stack trace corresponds to a plurality of prestack source and receiver location pairs. The prestack source and receiver location pairs may be associated with common depth point (CDP) traces which are used to create the stack trace. However, it is noted that the prestack source and receiver pairs are not necessarily associated with CDP traces. The plurality of stack traces may be organized into a 3-D stacked volume, referred to as a xe2x80x9ccubexe2x80x9d, as is well known in the art. Alternatively, the stacked data may be received, i.e., the received seismic data may already be stacked.
One important feature of numerous embodiments is that 3-D areal prestack data are generated from the stack traces. In a preferred embodiment, prestack gathers may be generated from each stacked trace by inverse moveout according to an azimuthally-dependent stacking velocity field, e.g., using a hyperbolic fanout of the stacking velocities, ray tracing seismic signals with a velocity and reflector model, or any other means of performing inverse moveout of the stacked trace.
Then, Delft 3-D multiple prediction may be applied to the generated 3-D areal prestack data to generate predicted multiples, as is well known in the art. The principle of reciprocity may be used to interchange the source and receiver locations of these generated 3-D areal prestack data whenever it may be convenient for data or computational management.
Finally, the results of the above process, i.e., the predicted multiples, may be output. For example, the predicted multiples may be output to storage, to an external system, e.g., to another computer system over a network, or to a display device, such as a computer monitor or printer.
It should be noted that simply applying the Delft technique in 3-D may be computationally expensive. For example, in 2-D, the Delft computations can require days to complete on a single computer system, and may require hours to compute on a large PC cluster. The 3-D computations could easily be a million times more compute-intensive than the 2-D computations. Thus, a number of tradeoffs or compromises may be made to allow useful 3-D multiple predictions to be made in a small fraction of the time otherwise required.
As an example, in one embodiment, each stacked trace may be fanned out only to a user-specified maximum offset. Said another way, the user may specify a maximum fanout distance from the location of each stacked trace, thereby specifying an aperture or neighborhood around the location to be used in the process. Since the computational cost is proportional to the square of this distance, significant time savings may be attained.
Another example of a compromise or restriction is to compute only the zero-offset traces in the 3-D Delft multiple prediction step, instead of computing all offsets and azimuths and then stacking the resulting prestack multiple predictions, thereby substantially reducing the number of required computations.
Although each of the above restrictions generally reduces one""s ability to accurately predict some classes of free-surface multiples, the tradeoffs are justifiable when the key multiples of interest are surface peglegs or pure surface multiples generated from moderately-dipping subsurface reflectors.
Under the restrictions described above, the algorithm for 3-D prestack/poststack multiple prediction is straightforward. 3D stacked seismic data may be received, and, at each common-depth-point (CDP) location, inverse moveout performed on the stacked trace (at that location) to all locations within a user-defined aperture, i.e., an area within a user-specified distance. In various embodiments, the specified area may have different geometric shapes, e.g., circle, ellipse, square, hexagon, etc., although for programming convenience and ease, a square or box is preferably used. It should be noted that the term xe2x80x9ccommon-depth-pointxe2x80x9d is often used in the art to denote xe2x80x9ccommon midpointxe2x80x9d or CMP, which technically refers to a CDP where the subsurface is horizontal. As used herein, the terms CDP and CMP may be used interchangeably. This shortcut also provides sizable computational savings over the theoretical cost of a full 3-D prestack Delft computation.
Inverse moveout may be performed in any of a variety of ways. In a preferred embodiment, inverse moveout may comprise inverse normal moveout (inverse NMO), as is well known in the art, where each stacked trace is modified (compressed and stretched) and translated to offset locations around the stacked trace location using shifted stacking hyperbolas, i.e., using stacking velocities. Other, more sophisticated, and computationally intensive, inverse moveout techniques contemplated include using a reflector and velocity model to ray trace arrivals of reflectors for inverse moveout parameterization, among others. In an alternate embodiment, stacking velocity and reflector times may be used to synthesize a spike trace at the same location as the selected stack trace, and inverse moveout applied to the selected stack trace and the spike trace to generate respective inverse moveout corrected traces at the stack seismic trace location.
Then, each inverse moveout corrected trace is convolved. For example, in one embodiment, autoconvolution or retro-correlation may be applied to the trace. In other words, the inverse moveout corrected trace may be convolved with itself. In an embodiment where a spike trace is synthesized and inverse moveout applied, the two inverse moveout corrected traces may be convolved with one another to generate a convolved trace.
The convolved traces are then summed at each CDP location. In the small aperture limit, 3-D prestack/poststack multiple prediction reduces to taking each stack trace and convolving it, either with itself or with a synthesized trace, in order to predict multiples. It is noted that autoconvolution of each stack trace is the classic 1966 approach of Anstey and Newman for 1-D prediction of multiples, and thus 3-D prestack/poststack multiple prediction may be considered to be an extension of their approach to 3-D. For further information regarding the approach of Anstey and Newman for 1-D prediction of multiples, please see xe2x80x9cThe Sectional Autocorrelogram and the Sectional Retrocorrelogramxe2x80x9d by Anstey and Newman, Geophysical Prospecting, v. 14, pp. 389-426, which is hereby incorporated by reference in its entirety, as though fully and completely set forth herein.
In another approach, prestack predicted multiple gathers may be created, for example, by generating a predicted multiple trace for each source-receiver pair of interest in the prestack data geometry. In one embodiment, the predicted multiple trace for each source-receiver pair may be generated as follows:
The predicted multiple trace may be initialized, e.g., with zero values, and an aperture selected, typically containing a neighborhood around both the source and receiver locations. For example, in one embodiment, the aperture or neighborhood may be an ellipse with the source at one focus and the receiver at another, although other geometries are also contemplated.
For each location in that aperture, inverse moveout may be applied to a stack trace from the source location to that location, and to a stack trace from the receiver location to that location. The two inverse moved out traces may then be convolved together and summed into the predicted multiple trace. The generated predicted multiple traces together form or compose the prestack predicted multiple gathers.
In variants of the above approach, an inverse moved out spike trace may be used instead of the inverse moved out stack trace at the receiver or source location.
It should be noted that in various embodiments, the individual convolved traces may optionally be scaled and/or filtered before summing. Such xe2x80x9cfold compensationxe2x80x9d, xe2x80x9crho-filterxe2x80x9d and xe2x80x9cobliquityxe2x80x9d adjustments are well known to those skilled in the relevant art. In a preferred embodiment, for speed considerations, subsequent adaptive subtraction (or visual overlay) may be relied upon to accommodate distortions that may arise when the more expensive adjustments are not applied.
Thus, based on stacked seismic data, estimates may be computed of prestack data that would have been recorded at any desired source/receiver/azimuth/offset (distance).
The generated predicted multiples may then be used in a variety of ways. For example, in one embodiment, a seismic data set and predicted multiples may be received. The seismic data set may be any seismic data that has the same source/receiver geometry as the predicted multiples. Typically, the seismic data have been acquired over a formation of interest, and preferably comprise unmigrated seismic data. In one embodiment, the seismic data set may be the received stacked seismic data described above, i.e., the 3D stacked volume. In another embodiment, the seismic data set may be prestack seismic data, e.g., original seismic data from which the 3D stacked volume was generated. In yet another embodiment, the seismic data set may be DMO-corrected data.
A transform may optionally be applied to the seismic data set and predicted multiples, thereby generating a transformed seismic data set and transformed predicted multiples, where the transformed predicted multiples and the transformed seismic data set are in a format suitable for comparison. The transform, applied to both the seismic data set and the predicted multiples, may be any type of transform desired. For example, if the seismic data set is prestack data and the predicted multiples are predicted prestack traces, the transform may be a stacking process, resulting in stacked seismic data and stacked multiples. As another example, the transform may be a migration process, as is well known in the art. Thus, the transform may comprise any operation or process on the seismic data and multiples that produces a useful result, some examples of which are described below.
In one embodiment, the (optionally transformed) predicted multiples may be adaptively subtracted from the (optionally transformed) seismic data set, thereby generating a processed seismic data set, which may then be output. The processed seismic data set thus preferably comprises the original seismic data set with the multiples substantially removed or filtered out. Thus, the processed seismic data set may more clearly indicate characteristics of the formation or subsurface structure of interest.
In another embodiment, the (optionally transformed) predicted multiples may be overlaid on the (optionally transformed) seismic data set to generate a multiples template, which may then be output. As is well known in the art, the multiples template may be used in a variety of ways to assist analysis of the formation or subsurface structure of interest. For example, the multiples template may indicate which seismic features in the seismic data are likely due to multiple reflections and thereby not true primary reflections from the formation of interest.
Thus, various embodiments may use various sources and forms of seismic data and velocity information, including field data and models, to predict 3-D multiples. These predicted multiples may then be used to analyze or further process seismic data, e.g., to discover or characterize subsurface formations.